%本函数用于在原有构型config0的基础上在x轴方向上添加一个原子。
%添加方式为：在不考虑新原子与其他原子的距离均大于原子直径的条件下，先确定原有构型的空间边界，再将这个空间离散化为100*100*100网格，将原子放入这些网格中心，取体系能量最小
function config_min = particle_addition(config0)
    dis_num = 100; %离散化为dis_num^3个网格
    N0 = length(config0);
    N = N0 + 1;
    x0 = zeros(1,N0);
    y0 = zeros(1,N0);
    z0 = zeros(1,N0);
    
    config = cell(1,N0+1);
    for i=1:N0
        config{i} = config0{i};
    end
    
    %en_min记录最小能量对应的构型{1}和能量{2}
    en_min = cell(1,2);
    en_min{1} = config0;
    en_min{2} = 0;

    for j=1:N0
        x0(j) = config0{j}(1);
        y0(j) = config0{j}(2);
        z0(j) = config0{j}(3);
    end
    
    %确定空间边界
    ext = 1.5; %在原有边界上进一步拓宽
    x_min = min(x0) - ext; x_max = max(x0) + ext;
    x_gap = (x_max - x_min)/dis_num;
    y_min = min(y0) - ext; y_max = max(y0) + ext;
    y_gap = (y_max - y_min)/dis_num;
    z_min = min(z0) - ext; z_max = max(z0) + ext;
    z_gap = (z_max - z_min)/dis_num;
    
    xlist = linspace(x_min + x_gap/2, x_max - x_gap/2, dis_num);
    ylist = linspace(y_min + y_gap/2, y_max - y_gap/2, dis_num);
    zlist = linspace(z_min + z_gap/2, z_max - z_gap/2, dis_num);
    
    for x=xlist
        for y=ylist
            for z=zlist
                pos = [x,y,z];
                config{N} = pos;
                en = sim_en_cal(config);
                if(en < en_min{2})
                    en_min{1} = config;
                    en_min{2} = en;
                end
            end
        end
    end
    config_min = en_min{1};
    return
end

%计算新原子与其他原子的相互作用势之和
function en = sim_en_cal(config)
    N = length(config);
    en = 0;
    for i=1:(N-1)
        en = en + LJ_pot(config{i},config{N});
    end
    return
end

%取原子直径为0.88sigma，截断半径为3.0sigma
%pt1,pt2均为1*3的数组，记录原子中心在三维直角坐标系下的3个坐标分量
function pot = LJ_pot(pt1, pt2)
    %根据pt1,pt2两点的坐标计算两点距离
    rvector = pt1 - pt2;
    r_square = sum(rvector.^2);
    r = sqrt(r_square);
    pot = r^(-12)-r^(-6); 
    
%     if(r < 0.88)
%         pot = 0.88^(-12)-0.88^(-6);
%     elseif(r > 3)
%         pot = 0;
%     else
%         pot = r^(-12)-r^(-6);   
%     end

    return
end

%旧版本：直接沿x轴添加原子，使得新原子不与其他原子发生体积重叠
% d = 0.88;
% N_old = length(config0);
% config = cell(1,N_old+1);
% 
% for i=1:N_old
%     config{i} = config0{i};
% end
% 
% new_pos = [0,0,0];
% 
% while 1
%     new_pos = new_pos + [d,0,0];
%     for j=1:N_old
%         rj_vector = config0{j} - new_pos;
%         rj_square = sum(rj_vector.^2);
%         rj = sqrt(rj_square);
%         if(rj < d)
%             judge = 0;
%             break
%         else
%             judge = 1;
%         end
%     end
%     if(judge == 1)
%         break
%     else
%     end
% end
% config{N_old+1} = new_pos;